<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  
  <title>OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之 kineticTheoryModel | Giskard&#39;s CFD Learning Tricks</title>
  <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1">
  <meta name="description" content="OpenFOAM 中双流体模型的 kineticTheoryModel 是以 “ Derivation, implementation, and validation of computer simulation models for gas-solid fluidized beds, B.G.M. van Wachem, Ph.D. Thesis, Delft University of Tec">
<meta property="og:type" content="article">
<meta property="og:title" content="OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之 kineticTheoryModel">
<meta property="og:url" content="http://xiaopingqiu.github.io/2015/09/19/kineticTheoryModel/index.html">
<meta property="og:site_name" content="Giskard's CFD Learning Tricks">
<meta property="og:description" content="OpenFOAM 中双流体模型的 kineticTheoryModel 是以 “ Derivation, implementation, and validation of computer simulation models for gas-solid fluidized beds, B.G.M. van Wachem, Ph.D. Thesis, Delft University of Tec">
<meta name="twitter:card" content="summary">
<meta name="twitter:title" content="OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之 kineticTheoryModel">
<meta name="twitter:description" content="OpenFOAM 中双流体模型的 kineticTheoryModel 是以 “ Derivation, implementation, and validation of computer simulation models for gas-solid fluidized beds, B.G.M. van Wachem, Ph.D. Thesis, Delft University of Tec">
  
  
    <link rel="icon" href="/favicon.png">
  
  <link href="//fonts.googleapis.com/css?family=Source+Code+Pro" rel="stylesheet" type="text/css">
  <link rel="stylesheet" href="/css/style.css" type="text/css">
  
<!-- Google Analytics -->
<script type="text/javascript">
(function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
(i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
})(window,document,'script','//www.google-analytics.com/analytics.js','ga');

##ga('create', '[object Object]', 'auto');
ga('create', 'UA-62501686-1', 'auto');
ga('send', 'pageview');

</script>
<!-- End Google Analytics -->


</head>
<body>
  <div id="container">
    <div id="wrap">
      <header id="header">
  <div id="banner"></div>
  <div id="header-outer" class="outer">
    <div id="header-title" class="inner">
      <h1 id="logo-wrap">
        <a href="/" id="logo">Giskard&#39;s CFD Learning Tricks</a>
      </h1>
      
        <h2 id="subtitle-wrap">
          <a href="/" id="subtitle">CFD and Scientific Computing</a>
        </h2>
      
    </div>
    <div id="header-inner" class="inner">
      <nav id="main-nav">
        <a id="main-nav-toggle" class="nav-icon"></a>
        
          <a class="main-nav-link" href="/">Home</a>
        
          <a class="main-nav-link" href="/archives">Archives</a>
        
          <a class="main-nav-link" href="/atom.xml">Rss</a>
        
          <a class="main-nav-link" href="/about">About</a>
        
      </nav>
      <nav id="sub-nav">
        
        <a id="nav-search-btn" class="nav-icon" title="Search"></a>
      </nav>
      <div id="search-form-wrap">
        <form action="//google.com/search" method="get" accept-charset="UTF-8" class="search-form"><input type="search" name="q" results="0" class="search-form-input" placeholder="Search"><button type="submit" class="search-form-submit">&#xF002;</button><input type="hidden" name="q" value="site:http://xiaopingqiu.github.io"></form>
      </div>
    </div>
  </div>
</header>

      <div class="outer">
        <section id="main"><article id="post-kineticTheoryModel" class="article article-type-post" itemscope itemprop="blogPost">
  <div class="article-meta">
    <a href="/2015/09/19/kineticTheoryModel/" class="article-date">
  <time datetime="2015-09-19T02:26:09.000Z" itemprop="datePublished">2015-09-19</time>
</a>
    
  <div class="article-category">
    <a class="article-category-link" href="/categories/OpenFOAM/">OpenFOAM</a>
  </div>

  </div>
  <div class="article-inner">
    
    
      <header class="article-header">
        
  
    <h1 class="article-title" itemprop="name">
      OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之 kineticTheoryModel
    </h1>
  

      </header>
    
    <div class="article-entry" itemprop="articleBody">
      
        <p>OpenFOAM 中双流体模型的 kineticTheoryModel 是以 “ <em>Derivation, implementation, and validation of computer simulation models for gas-solid fluidized beds</em>, B.G.M. van Wachem, Ph.D. Thesis, Delft University of Technology, Amsterdam, 2000. “ 为蓝本来设计的，下面分析这个类的代码。需要注意的是，这个类需要调用一些别的类（如 viscosityModel 等，后面会一一分析）来完成其功能。</p>
<a id="more"></a>
<h4 id="1-_头文件_kineticTheoryModel-H">1. 头文件 kineticTheoryModel.H</h4><p>头文件中要注意的是 kineticTheoryModel 类的继承关系，以及六个子类的智能指针作为 kineticTheoryModel 类的数据成员。</p>
<figure class="highlight lasso"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br><span class="line">28</span><br><span class="line">29</span><br><span class="line">30</span><br><span class="line">31</span><br><span class="line">32</span><br><span class="line">33</span><br><span class="line">34</span><br></pre></td><td class="code"><pre><span class="line">class kineticTheoryModel</span><br><span class="line">:</span><br><span class="line">    <span class="keyword">public</span> eddyViscosity</span><br><span class="line">    <span class="subst">&lt;</span></span><br><span class="line">        RASModel<span class="subst">&lt;</span>PhaseCompressibleTurbulenceModel<span class="subst">&lt;</span>phaseModel<span class="subst">&gt;</span> <span class="subst">&gt;</span></span><br><span class="line">    <span class="subst">&gt;</span> <span class="comment">// 继承自湍流类，所以这里 kineticTheoryModel 跟湍流类使用同样的接口</span></span><br><span class="line">&#123;</span><br><span class="line">    <span class="comment">// Private data</span></span><br><span class="line"></span><br><span class="line">        <span class="comment">// Input Fields</span></span><br><span class="line">            const phaseModel<span class="subst">&amp;</span> phase_;</span><br><span class="line">            </span><br><span class="line">        <span class="comment">// Sub-models </span></span><br><span class="line">	<span class="comment">// 下面五个 sub models 是kineticTheoryModel类为了实现其功能需要调用的类。</span></span><br><span class="line">	<span class="comment">//子类的智能指针定义为当前类的数据成员。注意这里的 "kineticTheoryModels" 是 namespace。</span></span><br><span class="line">        </span><br><span class="line">            <span class="comment">//- Run-time selected viscosity model</span></span><br><span class="line">            autoPtr<span class="subst">&lt;</span>kineticTheoryModels<span class="tag">::viscosityModel</span><span class="subst">&gt;</span> viscosityModel_;</span><br><span class="line">            </span><br><span class="line">            <span class="comment">//- Run-time selected conductivity model</span></span><br><span class="line">            autoPtr<span class="subst">&lt;</span>kineticTheoryModels<span class="tag">::conductivityModel</span><span class="subst">&gt;</span> conductivityModel_;</span><br><span class="line">            <span class="comment">//- Run-time selected radial distribution model</span></span><br><span class="line">            autoPtr<span class="subst">&lt;</span>kineticTheoryModels<span class="tag">::radialModel</span><span class="subst">&gt;</span> radialModel_;</span><br><span class="line"></span><br><span class="line">            <span class="comment">//- Run-time selected granular pressure model</span></span><br><span class="line">            autoPtr<span class="subst">&lt;</span>kineticTheoryModels<span class="tag">::granularPressureModel</span><span class="subst">&gt;</span></span><br><span class="line">                granularPressureModel_;</span><br><span class="line"></span><br><span class="line">            <span class="comment">//- Run-time selected frictional stress model</span></span><br><span class="line">            autoPtr<span class="subst">&lt;</span>kineticTheoryModels<span class="tag">::frictionalStressModel</span><span class="subst">&gt;</span></span><br><span class="line">                frictionalStressModel_;</span><br><span class="line"></span><br><span class="line">  <span class="attribute">...</span><span class="attribute">...</span></span><br><span class="line">  <span class="attribute">...</span><span class="attribute">...</span></span><br></pre></td></tr></table></figure>
<h4 id="2-_构造函数">2. 构造函数</h4><p>注意这里向基类传递的参数，这里的湍流类的继承关系比较复杂，比单相湍流复杂很多，后面会有具体的分析<br><figure class="highlight cpp"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br><span class="line">28</span><br><span class="line">29</span><br><span class="line">30</span><br><span class="line">31</span><br><span class="line">32</span><br><span class="line">33</span><br><span class="line">34</span><br><span class="line">35</span><br><span class="line">36</span><br><span class="line">37</span><br><span class="line">38</span><br><span class="line">39</span><br><span class="line">40</span><br><span class="line">41</span><br><span class="line">42</span><br><span class="line">43</span><br><span class="line">44</span><br><span class="line">45</span><br><span class="line">46</span><br><span class="line">47</span><br><span class="line">48</span><br><span class="line">49</span><br><span class="line">50</span><br><span class="line">51</span><br><span class="line">52</span><br><span class="line">53</span><br><span class="line">54</span><br><span class="line">55</span><br><span class="line">56</span><br><span class="line">57</span><br><span class="line">58</span><br><span class="line">59</span><br><span class="line">60</span><br><span class="line">61</span><br><span class="line">62</span><br><span class="line">63</span><br><span class="line">64</span><br></pre></td><td class="code"><pre><span class="line">Foam::RASModels::kineticTheoryModel::kineticTheoryModel</span><br><span class="line">(</span><br><span class="line">    <span class="keyword">const</span> volScalarField&amp; alpha,</span><br><span class="line">    <span class="keyword">const</span> volScalarField&amp; rho,</span><br><span class="line">    <span class="keyword">const</span> volVectorField&amp; U,</span><br><span class="line">    <span class="keyword">const</span> surfaceScalarField&amp; alphaRhoPhi,</span><br><span class="line">    <span class="keyword">const</span> surfaceScalarField&amp; phi,</span><br><span class="line">    <span class="keyword">const</span> transportModel&amp; phase,</span><br><span class="line">    <span class="keyword">const</span> word&amp; propertiesName,</span><br><span class="line">    <span class="keyword">const</span> word&amp; type</span><br><span class="line">)</span><br><span class="line">:</span><br><span class="line">    eddyViscosity&lt;RASModel&lt;PhaseCompressibleTurbulenceModel&lt;phaseModel&gt; &gt; &gt;</span><br><span class="line">    (</span><br><span class="line">        type,</span><br><span class="line">        alpha,</span><br><span class="line">        rho,</span><br><span class="line">        U,</span><br><span class="line">        alphaRhoPhi,</span><br><span class="line">        phi,</span><br><span class="line">        phase,</span><br><span class="line">        propertiesName</span><br><span class="line">    ),</span><br><span class="line"></span><br><span class="line">    phase_(phase),</span><br><span class="line"></span><br><span class="line">    viscosityModel_</span><br><span class="line">    (</span><br><span class="line">        kineticTheoryModels::viscosityModel::New</span><br><span class="line">        (</span><br><span class="line">            <span class="keyword">this</span>-&gt;coeffDict_ <span class="comment">// coeffDict_ 是 RASModel的成员</span></span><br><span class="line">        )</span><br><span class="line">    ),</span><br><span class="line">    conductivityModel_</span><br><span class="line">    (</span><br><span class="line">        kineticTheoryModels::conductivityModel::New</span><br><span class="line">        (</span><br><span class="line">            <span class="keyword">this</span>-&gt;coeffDict_</span><br><span class="line">        )</span><br><span class="line">    ),</span><br><span class="line">    radialModel_</span><br><span class="line">    (</span><br><span class="line">        kineticTheoryModels::radialModel::New</span><br><span class="line">        (</span><br><span class="line">            <span class="keyword">this</span>-&gt;coeffDict_</span><br><span class="line">        )</span><br><span class="line">    ),</span><br><span class="line">    granularPressureModel_</span><br><span class="line">    (</span><br><span class="line">        kineticTheoryModels::granularPressureModel::New</span><br><span class="line">        (</span><br><span class="line">            <span class="keyword">this</span>-&gt;coeffDict_</span><br><span class="line">        )</span><br><span class="line">    ),</span><br><span class="line">    frictionalStressModel_</span><br><span class="line">    (</span><br><span class="line">        kineticTheoryModels::frictionalStressModel::New</span><br><span class="line">        (</span><br><span class="line">            <span class="keyword">this</span>-&gt;coeffDict_</span><br><span class="line">        )</span><br><span class="line">    ),</span><br><span class="line"></span><br><span class="line">......</span><br><span class="line">......</span><br></pre></td></tr></table></figure></p>
<h4 id="3-_主要成员函数">3. 主要成员函数</h4><figure class="highlight gherkin"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br></pre></td><td class="code"><pre><span class="line">// <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> Member Functions  <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> //</span><br><span class="line"></span><br><span class="line">Foam::tmp<span class="variable">&lt;Foam::volScalarField&gt;</span></span><br><span class="line">Foam::RASModels::kineticTheoryModel::k() const</span><br><span class="line">&#123;</span><br><span class="line">    // notImplemented 是 Error.H中定义的一个函数，用于提示某个模型或函数没有实现。</span><br><span class="line">    notImplemented(<span class="string">"kineticTheoryModel::k()"</span>); </span><br><span class="line">    return nut_;</span><br><span class="line">&#125;</span><br><span class="line"></span><br><span class="line"></span><br><span class="line">Foam::tmp<span class="variable">&lt;Foam::volScalarField&gt;</span></span><br><span class="line">Foam::RASModels::kineticTheoryModel::epsilon() const</span><br><span class="line">&#123;</span><br><span class="line">    notImplemented(<span class="string">"kineticTheoryModel::epsilon()"</span>);</span><br><span class="line">    return nut_;</span><br><span class="line">&#125;</span><br></pre></td></tr></table></figure>
<p> <code>k</code> 和  <code>epsilon</code> 是两个从基类继承下来的函数，这两个函数在基类中是纯虚函数，所以虽然 kineticTheoryModel 用不到它们，但是还是需要进行定义，否则 kineticTheoryModel 类就将是一个虚基类，从而无法创建对象了。</p>
<figure class="highlight elixir"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br></pre></td><td class="code"><pre><span class="line"><span class="constant">Foam:</span><span class="symbol">:tmp&lt;Foam</span><span class="symbol">:</span><span class="symbol">:volSymmTensorField&gt;</span></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:RASModels</span><span class="symbol">:</span><span class="symbol">:kineticTheoryModel</span><span class="symbol">:</span><span class="symbol">:R</span>() const</span><br><span class="line">&#123;</span><br><span class="line">    <span class="keyword">return</span> tmp&lt;volSymmTensorField&gt;</span><br><span class="line">    (</span><br><span class="line">        new volSymmTensorField</span><br><span class="line">        (</span><br><span class="line">            <span class="constant">IOobject</span></span><br><span class="line">            (</span><br><span class="line">                <span class="constant">IOobject:</span><span class="symbol">:groupName</span>(<span class="string">"R"</span>, this-&gt;<span class="constant">U_.</span>group()),</span><br><span class="line">                this-&gt;runTime<span class="constant">_</span>.timeName(),</span><br><span class="line">                this-&gt;mesh<span class="constant">_</span>,</span><br><span class="line">                <span class="constant">IOobject:</span><span class="symbol">:NO_READ</span>,</span><br><span class="line">                <span class="constant">IOobject:</span><span class="symbol">:NO_WRITE</span></span><br><span class="line">            ),</span><br><span class="line">          - (this-&gt;nut<span class="constant">_</span>)*dev(twoSymm(<span class="symbol">fvc:</span><span class="symbol">:grad</span>(this-&gt;<span class="constant">U_)</span>))</span><br><span class="line">          - (lambda<span class="constant">_</span>*<span class="symbol">fvc:</span><span class="symbol">:div</span>(this-&gt;phi<span class="constant">_</span>))*<span class="symbol">symmTensor:</span><span class="symbol">:I</span></span><br><span class="line">        )</span><br><span class="line">    );</span><br><span class="line">&#125;</span><br></pre></td></tr></table></figure>
<p>R 函数返回所谓的雷诺应力，其实跟单相湍流类的 R 函数返回值有点差别：<br>$$<br>-\nu_t[\nabla U + \nabla U^T-\frac{2}{3}(\nabla \cdot U)\cdot \mathrm{I}] - \lambda(\nabla \cdot U)\cdot \mathrm{I}<br>$$<br>这里的 R 与下面的 devRoReff 函数返回值只相差一个 rho。<br>而单相湍流模型的 R 则为：<br>$$<br>\frac{2}{3}k\cdot \mathrm{I}- \nu_t (\nabla U + \nabla U^T)<br>$$</p>
<figure class="highlight elixir"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br><span class="line">28</span><br><span class="line">29</span><br><span class="line">30</span><br><span class="line">31</span><br><span class="line">32</span><br><span class="line">33</span><br><span class="line">34</span><br><span class="line">35</span><br><span class="line">36</span><br><span class="line">37</span><br><span class="line">38</span><br><span class="line">39</span><br><span class="line">40</span><br><span class="line">41</span><br><span class="line">42</span><br><span class="line">43</span><br><span class="line">44</span><br><span class="line">45</span><br><span class="line">46</span><br><span class="line">47</span><br><span class="line">48</span><br><span class="line">49</span><br><span class="line">50</span><br><span class="line">51</span><br><span class="line">52</span><br><span class="line">53</span><br><span class="line">54</span><br><span class="line">55</span><br></pre></td><td class="code"><pre><span class="line"></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:tmp&lt;Foam</span><span class="symbol">:</span><span class="symbol">:volScalarField&gt;</span></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:RASModels</span><span class="symbol">:</span><span class="symbol">:kineticTheoryModel</span><span class="symbol">:</span><span class="symbol">:pPrime</span>() const</span><br><span class="line">&#123;</span><br><span class="line">    <span class="regexp">//</span> <span class="constant">Local </span>references</span><br><span class="line">    const volScalarField&amp; alpha = this-&gt;alpha<span class="constant">_</span>;</span><br><span class="line">    const volScalarField&amp; rho = phase<span class="constant">_</span>.rho();</span><br><span class="line"></span><br><span class="line">    <span class="keyword">return</span></span><br><span class="line">    (</span><br><span class="line">        <span class="constant">Theta_</span></span><br><span class="line">       *granularPressureModel<span class="constant">_</span>-&gt;granularPressureCoeffPrime</span><br><span class="line">        (</span><br><span class="line">            alpha,</span><br><span class="line">            radialModel<span class="constant">_</span>-&gt;g<span class="number">0</span>(alpha, alphaMinFriction<span class="constant">_</span>, alphaMax<span class="constant">_</span>),</span><br><span class="line">            radialModel<span class="constant">_</span>-&gt;g0prime(alpha, alphaMinFriction<span class="constant">_</span>, alphaMax<span class="constant">_</span>),</span><br><span class="line">            rho,</span><br><span class="line">            e<span class="constant">_</span></span><br><span class="line">        )</span><br><span class="line">     +  frictionalStressModel<span class="constant">_</span>-&gt;frictionalPressurePrime</span><br><span class="line">        (</span><br><span class="line">            alpha,</span><br><span class="line">            alphaMinFriction<span class="constant">_</span>,</span><br><span class="line">            alphaMax<span class="constant">_</span></span><br><span class="line">        )</span><br><span class="line">    );</span><br><span class="line">&#125;</span><br><span class="line"></span><br><span class="line"></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:tmp&lt;Foam</span><span class="symbol">:</span><span class="symbol">:surfaceScalarField&gt;</span></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:RASModels</span><span class="symbol">:</span><span class="symbol">:kineticTheoryModel</span><span class="symbol">:</span><span class="symbol">:pPrimef</span>() const</span><br><span class="line">&#123;</span><br><span class="line">    <span class="regexp">//</span> <span class="constant">Local </span>references</span><br><span class="line">    const volScalarField&amp; alpha = this-&gt;alpha<span class="constant">_</span>;</span><br><span class="line">    const volScalarField&amp; rho = phase<span class="constant">_</span>.rho();</span><br><span class="line"></span><br><span class="line">    <span class="keyword">return</span> <span class="symbol">fvc:</span><span class="symbol">:interpolate</span></span><br><span class="line">    (</span><br><span class="line">        <span class="constant">Theta_</span></span><br><span class="line">       *granularPressureModel<span class="constant">_</span>-&gt;granularPressureCoeffPrime</span><br><span class="line">        (</span><br><span class="line">            alpha,</span><br><span class="line">            radialModel<span class="constant">_</span>-&gt;g<span class="number">0</span>(alpha, alphaMinFriction<span class="constant">_</span>, alphaMax<span class="constant">_</span>),</span><br><span class="line">            radialModel<span class="constant">_</span>-&gt;g0prime(alpha, alphaMinFriction<span class="constant">_</span>, alphaMax<span class="constant">_</span>),</span><br><span class="line">            rho,</span><br><span class="line">            e<span class="constant">_</span></span><br><span class="line">        )</span><br><span class="line">     +  frictionalStressModel<span class="constant">_</span>-&gt;frictionalPressurePrime</span><br><span class="line">        (</span><br><span class="line">            alpha,</span><br><span class="line">            alphaMinFriction<span class="constant">_</span>,</span><br><span class="line">            alphaMax<span class="constant">_</span></span><br><span class="line">        )</span><br><span class="line">    );</span><br><span class="line">&#125;</span><br></pre></td></tr></table></figure>
<p>pPrime 和 pPrimef 两个函数，返回的是固相压力对固相孔隙率的导数（$\partial P_s/\partial \varepsilon_s$）。</p>
<p>两个在”UEqn.H” 被动量方程调用的函数<br><figure class="highlight elixir"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br><span class="line">28</span><br><span class="line">29</span><br><span class="line">30</span><br><span class="line">31</span><br><span class="line">32</span><br><span class="line">33</span><br><span class="line">34</span><br><span class="line">35</span><br><span class="line">36</span><br><span class="line">37</span><br><span class="line">38</span><br><span class="line">39</span><br><span class="line">40</span><br></pre></td><td class="code"><pre><span class="line"><span class="constant">Foam:</span><span class="symbol">:tmp&lt;Foam</span><span class="symbol">:</span><span class="symbol">:volSymmTensorField&gt;</span></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:RASModels</span><span class="symbol">:</span><span class="symbol">:kineticTheoryModel</span><span class="symbol">:</span><span class="symbol">:devRhoReff</span>() const</span><br><span class="line">&#123;</span><br><span class="line">    <span class="keyword">return</span> tmp&lt;volSymmTensorField&gt;</span><br><span class="line">    (</span><br><span class="line">        new volSymmTensorField</span><br><span class="line">        (</span><br><span class="line">            <span class="constant">IOobject</span></span><br><span class="line">            (</span><br><span class="line">                <span class="constant">IOobject:</span><span class="symbol">:groupName</span>(<span class="string">"devRhoReff"</span>, this-&gt;<span class="constant">U_.</span>group()),</span><br><span class="line">                this-&gt;runTime<span class="constant">_</span>.timeName(),</span><br><span class="line">                this-&gt;mesh<span class="constant">_</span>,</span><br><span class="line">                <span class="constant">IOobject:</span><span class="symbol">:NO_READ</span>,</span><br><span class="line">                <span class="constant">IOobject:</span><span class="symbol">:NO_WRITE</span></span><br><span class="line">            ),</span><br><span class="line">          - (this-&gt;rho<span class="constant">_</span>*this-&gt;nut<span class="constant">_</span>)</span><br><span class="line">           *dev(twoSymm(<span class="symbol">fvc:</span><span class="symbol">:grad</span>(this-&gt;<span class="constant">U_)</span>))</span><br><span class="line">          - ((this-&gt;rho<span class="constant">_</span>*lambda<span class="constant">_</span>)*<span class="symbol">fvc:</span><span class="symbol">:div</span>(this-&gt;phi<span class="constant">_</span>))*<span class="symbol">symmTensor:</span><span class="symbol">:I</span></span><br><span class="line">        )</span><br><span class="line">    );</span><br><span class="line">&#125;</span><br><span class="line"></span><br><span class="line"></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:tmp&lt;Foam</span><span class="symbol">:</span><span class="symbol">:fvVectorMatrix&gt;</span></span><br><span class="line"><span class="constant">Foam:</span><span class="symbol">:RASModels</span><span class="symbol">:</span><span class="symbol">:kineticTheoryModel</span><span class="symbol">:</span><span class="symbol">:divDevRhoReff</span></span><br><span class="line">(</span><br><span class="line">    volVectorField&amp; <span class="constant">U</span></span><br><span class="line">) const</span><br><span class="line">&#123;</span><br><span class="line">    <span class="keyword">return</span></span><br><span class="line">    (</span><br><span class="line">      - <span class="symbol">fvm:</span><span class="symbol">:laplacian</span>(this-&gt;rho<span class="constant">_</span>*this-&gt;nut<span class="constant">_</span>, <span class="constant">U)</span></span><br><span class="line">      - <span class="symbol">fvc:</span><span class="symbol">:div</span></span><br><span class="line">        (</span><br><span class="line">            (this-&gt;rho<span class="constant">_</span>*this-&gt;nut<span class="constant">_</span>)*dev2(<span class="constant">T(</span><span class="symbol">fvc:</span><span class="symbol">:grad</span>(<span class="constant">U)</span>))</span><br><span class="line">          + ((this-&gt;rho<span class="constant">_</span>*lambda<span class="constant">_</span>)*<span class="symbol">fvc:</span><span class="symbol">:div</span>(this-&gt;phi<span class="constant">_</span>))</span><br><span class="line">           *dimensioned&lt;symmTensor&gt;(<span class="string">"I"</span>, dimless, <span class="symbol">symmTensor:</span><span class="symbol">:I</span>)</span><br><span class="line">        )</span><br><span class="line">    );</span><br><span class="line">&#125;</span><br></pre></td></tr></table></figure></p>
<p>对应的公式分别为：<br><strong>devRhoReff</strong><br>$$<br>- \rho  \nu_t \cdot dev(\nabla U + \nabla U^T) - \rho \lambda(\nabla \cdot U)\cdot  \mathrm{I} \\<br> = - \rho  \nu_t(\nabla U + \nabla U^T) +  \rho  \nu_t \frac{2}{3}(\nabla \cdot U)\cdot \mathrm{I} - \rho \lambda (\nabla \cdot U)\cdot  \mathrm{I}<br>$$</p>
<p><strong>divDevRhoReff</strong><br>$$<br>\ - \nabla \cdot (\rho \nu_t \nabla U) - \nabla \cdot \left [\rho \nu_t \nabla U^T\ - \rho \nu_t \frac{2}{3}(\nabla \cdot U)\cdot \mathrm{I} \right ] + \rho \lambda (\nabla \cdot U) \cdot \mathrm{I} )<br>$$</p>
<p> <code>correct</code>是计算颗粒温度 <code>Theta</code> 的函数，这是 kineticTheoryModel 类最重要的部分。<br><figure class="highlight stata"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br><span class="line">28</span><br><span class="line">29</span><br><span class="line">30</span><br><span class="line">31</span><br><span class="line">32</span><br><span class="line">33</span><br><span class="line">34</span><br><span class="line">35</span><br><span class="line">36</span><br><span class="line">37</span><br><span class="line">38</span><br><span class="line">39</span><br><span class="line">40</span><br><span class="line">41</span><br><span class="line">42</span><br><span class="line">43</span><br><span class="line">44</span><br><span class="line">45</span><br><span class="line">46</span><br><span class="line">47</span><br><span class="line">48</span><br><span class="line">49</span><br><span class="line">50</span><br><span class="line">51</span><br><span class="line">52</span><br><span class="line">53</span><br><span class="line">54</span><br><span class="line">55</span><br><span class="line">56</span><br><span class="line">57</span><br><span class="line">58</span><br><span class="line">59</span><br><span class="line">60</span><br><span class="line">61</span><br><span class="line">62</span><br><span class="line">63</span><br><span class="line">64</span><br><span class="line">65</span><br><span class="line">66</span><br><span class="line">67</span><br><span class="line">68</span><br><span class="line">69</span><br><span class="line">70</span><br><span class="line">71</span><br><span class="line">72</span><br><span class="line">73</span><br><span class="line">74</span><br><span class="line">75</span><br><span class="line">76</span><br><span class="line">77</span><br><span class="line">78</span><br><span class="line">79</span><br><span class="line">80</span><br><span class="line">81</span><br><span class="line">82</span><br><span class="line">83</span><br><span class="line">84</span><br><span class="line">85</span><br><span class="line">86</span><br><span class="line">87</span><br><span class="line">88</span><br><span class="line">89</span><br><span class="line">90</span><br><span class="line">91</span><br><span class="line">92</span><br><span class="line">93</span><br><span class="line">94</span><br><span class="line">95</span><br><span class="line">96</span><br><span class="line">97</span><br><span class="line">98</span><br><span class="line">99</span><br><span class="line">100</span><br><span class="line">101</span><br><span class="line">102</span><br><span class="line">103</span><br><span class="line">104</span><br><span class="line">105</span><br><span class="line">106</span><br><span class="line">107</span><br><span class="line">108</span><br><span class="line">109</span><br><span class="line">110</span><br><span class="line">111</span><br><span class="line">112</span><br><span class="line">113</span><br><span class="line">114</span><br><span class="line">115</span><br><span class="line">116</span><br><span class="line">117</span><br><span class="line">118</span><br><span class="line">119</span><br><span class="line">120</span><br><span class="line">121</span><br><span class="line">122</span><br><span class="line">123</span><br><span class="line">124</span><br><span class="line">125</span><br><span class="line">126</span><br><span class="line">127</span><br><span class="line">128</span><br><span class="line">129</span><br><span class="line">130</span><br><span class="line">131</span><br><span class="line">132</span><br><span class="line">133</span><br><span class="line">134</span><br><span class="line">135</span><br><span class="line">136</span><br><span class="line">137</span><br><span class="line">138</span><br><span class="line">139</span><br><span class="line">140</span><br><span class="line">141</span><br><span class="line">142</span><br><span class="line">143</span><br><span class="line">144</span><br><span class="line">145</span><br><span class="line">146</span><br><span class="line">147</span><br><span class="line">148</span><br><span class="line">149</span><br><span class="line">150</span><br><span class="line">151</span><br><span class="line">152</span><br><span class="line">153</span><br><span class="line">154</span><br><span class="line">155</span><br><span class="line">156</span><br><span class="line">157</span><br><span class="line">158</span><br><span class="line">159</span><br><span class="line">160</span><br><span class="line">161</span><br><span class="line">162</span><br><span class="line">163</span><br><span class="line">164</span><br><span class="line">165</span><br><span class="line">166</span><br><span class="line">167</span><br><span class="line">168</span><br><span class="line">169</span><br><span class="line">170</span><br><span class="line">171</span><br><span class="line">172</span><br><span class="line">173</span><br><span class="line">174</span><br><span class="line">175</span><br><span class="line">176</span><br><span class="line">177</span><br><span class="line">178</span><br><span class="line">179</span><br><span class="line">180</span><br><span class="line">181</span><br><span class="line">182</span><br><span class="line">183</span><br><span class="line">184</span><br><span class="line">185</span><br><span class="line">186</span><br><span class="line">187</span><br><span class="line">188</span><br><span class="line">189</span><br><span class="line">190</span><br><span class="line">191</span><br><span class="line">192</span><br><span class="line">193</span><br><span class="line">194</span><br><span class="line">195</span><br><span class="line">196</span><br><span class="line">197</span><br><span class="line">198</span><br><span class="line">199</span><br><span class="line">200</span><br><span class="line">201</span><br><span class="line">202</span><br><span class="line">203</span><br><span class="line">204</span><br></pre></td><td class="code"><pre><span class="line">void Foam::RASModels::kineticTheoryModel::correct()</span><br><span class="line">&#123;</span><br><span class="line">    <span class="comment">// Local references</span></span><br><span class="line">    volScalarField <span class="keyword">alpha</span>(<span class="literal">max</span>(this-&gt;alpha_, <span class="literal">scalar</span>(0)));</span><br><span class="line">    <span class="keyword">const</span> volScalarField&amp; rho = phase_.rho();</span><br><span class="line">    <span class="keyword">const</span> surfaceScalarField&amp; alphaRhoPhi = this-&gt;alphaRhoPhi_;</span><br><span class="line">    <span class="keyword">const</span> volVectorField&amp; <span class="keyword">U</span> = this-&gt;U_; <span class="comment">// 当前相的速度</span></span><br><span class="line">    <span class="keyword">const</span> volVectorField&amp; Uc_ = phase_.fluid().otherPhase(phase_).<span class="keyword">U</span>(); <span class="comment">// 另一相的速度</span></span><br><span class="line"></span><br><span class="line">    <span class="keyword">const</span> <span class="keyword">scalar</span> sqrtPi = <span class="literal">sqrt</span>(constant::mathematical::pi);</span><br><span class="line">    dimensionedScalar ThetaSmall(<span class="string">"ThetaSmall"</span>, Theta_.dimensions(), 1.0e-6);</span><br><span class="line">    dimensionedScalar ThetaSmallSqrt(<span class="literal">sqrt</span>(ThetaSmall));</span><br><span class="line"></span><br><span class="line">    tmp&lt;volScalarField&gt; tda(phase_.<span class="literal">d</span>()); <span class="comment">// 颗粒直径</span></span><br><span class="line">    <span class="keyword">const</span> volScalarField&amp; da = tda();</span><br><span class="line"></span><br><span class="line">    tmp&lt;volTensorField&gt; tgradU(fvc::grad(this-&gt;U_));</span><br><span class="line">    <span class="keyword">const</span> volTensorField&amp; gradU(tgradU());</span><br><span class="line">    volSymmTensorField <span class="literal">D</span>(symm(gradU));</span><br><span class="line"></span><br><span class="line">    <span class="comment">// Calculating the radial distribution function </span></span><br><span class="line">    <span class="comment">// 调用 radialModel 类 来计算径向分布函数</span></span><br><span class="line">    gs0_ = radialModel_-&gt;g0(<span class="keyword">alpha</span>, alphaMinFriction_, alphaMax_);</span><br><span class="line"></span><br><span class="line">    <span class="keyword">if</span> (!equilibrium_) <span class="comment">// 如果 equilibrium_ = off ，那么将用这个偏微分方程来计算颗粒温度，否则将改用一个代数方程，见下文</span></span><br><span class="line">    &#123;</span><br><span class="line">        <span class="comment">// Particle viscosity (Table 3.2, p.47)</span></span><br><span class="line">	<span class="comment">// 调用viscosityModel 来更新颗粒相的粘度</span></span><br><span class="line">        nut_ = viscosityModel_-&gt;nu(<span class="keyword">alpha</span>, Theta_, gs0_, rho, da, e_);</span><br><span class="line"></span><br><span class="line">        volScalarField ThetaSqrt(<span class="literal">sqrt</span>(Theta_));</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Bulk viscosity  p. 45 (Lun et al. 1984). </span></span><br><span class="line">        lambda_ = (4.0/3.0)*sqr(<span class="keyword">alpha</span>)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Stress tensor, Definitions, Table 3.1, p. 43</span></span><br><span class="line">        volSymmTensorField tau <span class="comment">// 颗粒应力</span></span><br><span class="line">        (</span><br><span class="line">            rho*(2.0*nut_*<span class="keyword">D</span> + (lambda_ - (2.0/3.0)*nut_)*tr(<span class="keyword">D</span>)*I)</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Dissipation (Eq. 3.24, p.50)</span></span><br><span class="line">        volScalarField gammaCoeff <span class="comment">//颗粒动能耗散项</span></span><br><span class="line">        (</span><br><span class="line">            12.0*(1.0 - sqr(e_))</span><br><span class="line"><span class="comment">           *max(sqr(alpha), residualAlpha_)</span></span><br><span class="line"><span class="comment">           *rho*gs0_*(1.0/da)*ThetaSqrt/sqrtPi</span></span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Drag // 调用曳力模型计算曳力系数</span></span><br><span class="line">        volScalarField beta(phase_.fluid().drag(phase_).K());</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Eq. 3.25, p. 50 Js = J1 - J2</span></span><br><span class="line">        volScalarField J1(3.0*beta);</span><br><span class="line">        volScalarField J2</span><br><span class="line">        (</span><br><span class="line">            0.25*sqr(beta)*da*magSqr(<span class="keyword">U</span> - Uc_)</span><br><span class="line">           /(</span><br><span class="line">               <span class="literal">max</span>(<span class="keyword">alpha</span>, residualAlpha_)*rho</span><br><span class="line"><span class="comment">              *sqrtPi*(ThetaSqrt + ThetaSmallSqrt)</span></span><br><span class="line">            )</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        <span class="comment">// particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)</span></span><br><span class="line">        volScalarField PsCoeff</span><br><span class="line">        (</span><br><span class="line">            granularPressureModel_-&gt;granularPressureCoeff</span><br><span class="line">            (</span><br><span class="line">                <span class="keyword">alpha</span>,</span><br><span class="line">                gs0_,</span><br><span class="line">                rho,</span><br><span class="line">                e_</span><br><span class="line">            )</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        <span class="comment">// 'thermal' conductivity (Table 3.3, p. 49)</span></span><br><span class="line">	<span class="comment">// 调用 conductivityModel 计算颗粒脉动能量的传导系数</span></span><br><span class="line">        kappa_ = conductivityModel_-&gt;<span class="keyword">kappa</span>(<span class="keyword">alpha</span>, Theta_, gs0_, rho, da, e_);</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Construct the granular temperature equation (Eq. 3.20, p. 44)</span></span><br><span class="line">        <span class="comment">// NB. note that there are two typos in Eq. 3.20:</span></span><br><span class="line">        <span class="comment">//     Ps should be without grad</span></span><br><span class="line">        <span class="comment">//     the laplacian has the wrong sign</span></span><br><span class="line">	<span class="comment">// 构建颗粒温度方程，注意，开头提到的文献里的颗粒温度方程有两处 typo，</span></span><br><span class="line">	<span class="comment">//下面的代码修复了这两处错误，后文会给出正确的公式。</span></span><br><span class="line">        fvScalarMatrix ThetaEqn</span><br><span class="line">        (</span><br><span class="line">            1.5*</span><br><span class="line">            (</span><br><span class="line">                fvm::ddt(<span class="keyword">alpha</span>, rho, Theta_)</span><br><span class="line">              + fvm::div(alphaRhoPhi, Theta_)</span><br><span class="line">              - fvc::Sp(fvc::ddt(<span class="keyword">alpha</span>, rho) + fvc::div(alphaRhoPhi), Theta_)</span><br><span class="line">            )</span><br><span class="line">          - fvm::laplacian(kappa_, Theta_, <span class="string">"laplacian(kappa,Theta)"</span>)</span><br><span class="line">         ==</span><br><span class="line">            fvm::SuSp(-((PsCoeff*I) &amp;&amp; gradU), Theta_)</span><br><span class="line">          + (tau &amp;&amp; gradU)</span><br><span class="line">          + fvm::Sp(-gammaCoeff, Theta_)</span><br><span class="line">          + fvm::Sp(-J1, Theta_)</span><br><span class="line">          + fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        ThetaEqn.relax();</span><br><span class="line">        ThetaEqn.solve();</span><br><span class="line">    &#125;</span><br><span class="line">    <span class="keyword">else</span> <span class="comment">// 如果 equilibrium = on， 将使用一个代数方程来计算颗粒温度。</span></span><br><span class="line">    &#123;</span><br><span class="line">        <span class="comment">// Equilibrium =&gt; dissipation == production</span></span><br><span class="line">        <span class="comment">// Eq. 4.14, p.82</span></span><br><span class="line">        volScalarField K1(2.0*(1.0 + e_)*rho*gs0_);</span><br><span class="line">        volScalarField K3</span><br><span class="line">        (</span><br><span class="line">            0.5*da*rho*</span><br><span class="line">            (</span><br><span class="line">                (sqrtPi/(3.0*(3.0 - e_)))</span><br><span class="line"><span class="comment">               *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_)</span></span><br><span class="line">               +1.6*<span class="keyword">alpha</span>*gs0_*(1.0 + e_)/sqrtPi</span><br><span class="line">            )</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        volScalarField K2</span><br><span class="line">        (</span><br><span class="line">            4.0*da*rho*(1.0 + e_)*<span class="keyword">alpha</span>*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        volScalarField K4(12.0*(1.0 - sqr(e_))*rho*gs0_/(da*sqrtPi));</span><br><span class="line"></span><br><span class="line">        volScalarField trD</span><br><span class="line">        (</span><br><span class="line">            <span class="keyword">alpha</span>/(<span class="keyword">alpha</span> + residualAlpha_)</span><br><span class="line"><span class="comment">           *fvc::div(this-&gt;phi_)</span></span><br><span class="line">        );</span><br><span class="line">        volScalarField tr2D(sqr(trD));</span><br><span class="line">        volScalarField trD2(tr(<span class="keyword">D</span> &amp; <span class="keyword">D</span>));</span><br><span class="line"></span><br><span class="line">        volScalarField t1(K1*<span class="keyword">alpha</span> + rho);</span><br><span class="line">        volScalarField l1(-t1*trD);</span><br><span class="line">        volScalarField l2(sqr(t1)*tr2D);</span><br><span class="line">        volScalarField l3</span><br><span class="line">        (</span><br><span class="line">            4.0</span><br><span class="line"><span class="comment">           *K4</span></span><br><span class="line"><span class="comment">           *alpha</span></span><br><span class="line"><span class="comment">           *(2.0*K3*trD2 + K2*tr2D)</span></span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        Theta_ = sqr</span><br><span class="line">        (</span><br><span class="line">            (l1 + <span class="literal">sqrt</span>(l2 + l3))</span><br><span class="line">           /(2.0*<span class="literal">max</span>(<span class="keyword">alpha</span>, residualAlpha_)*K4)</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        kappa_ = conductivityModel_-&gt;<span class="keyword">kappa</span>(<span class="keyword">alpha</span>, Theta_, gs0_, rho, da, e_);</span><br><span class="line">    &#125;</span><br><span class="line"></span><br><span class="line"></span><br><span class="line">   <span class="comment">// 限定 颗粒温度的上下限。</span></span><br><span class="line">  <span class="comment">// max 和 min 函数的定义，没有找到。经验证， max 的作用是让小于0的归零， min 是让大于100的等于100。</span></span><br><span class="line">    Theta_.<span class="literal">max</span>(0);</span><br><span class="line">    Theta_.<span class="literal">min</span>(100);</span><br><span class="line"></span><br><span class="line">    &#123;</span><br><span class="line">	<span class="comment">//利用先得到的颗粒温度更新颗粒相的粘度</span></span><br><span class="line">        <span class="comment">// particle viscosity (Table 3.2, p.47)</span></span><br><span class="line">        nut_ = viscosityModel_-&gt;nu(<span class="keyword">alpha</span>, Theta_, gs0_, rho, da, e_);</span><br><span class="line"></span><br><span class="line">        volScalarField ThetaSqrt(<span class="literal">sqrt</span>(Theta_));</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Bulk viscosity  p. 45 (Lun et al. 1984).</span></span><br><span class="line">        lambda_ = (4.0/3.0)*sqr(<span class="keyword">alpha</span>)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Frictional pressure // 计算由于颗粒之间的摩擦作用产生的一个等效的颗粒相压力作用。</span></span><br><span class="line">        volScalarField pf</span><br><span class="line">        (</span><br><span class="line">            frictionalStressModel_-&gt;frictionalPressure</span><br><span class="line">            (</span><br><span class="line">                <span class="keyword">alpha</span>,</span><br><span class="line">                alphaMinFriction_,</span><br><span class="line">                alphaMax_</span><br><span class="line">            )</span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Add frictional shear viscosity, Eq. 3.30, p. 52</span></span><br><span class="line">       <span class="comment">// 将颗粒摩擦产生的颗粒相粘度加到由颗粒温度计算得到的颗粒相粘度中，作为总的颗粒相粘度</span></span><br><span class="line">        nut_ += frictionalStressModel_-&gt;nu </span><br><span class="line">        (</span><br><span class="line">            <span class="keyword">alpha</span>,</span><br><span class="line">            alphaMax_,</span><br><span class="line">            pf/rho,</span><br><span class="line">            <span class="literal">D</span></span><br><span class="line">        );</span><br><span class="line"></span><br><span class="line">        <span class="comment">// Limit viscosity</span></span><br><span class="line">	<span class="comment">// 限定颗粒相粘度的上限</span></span><br><span class="line">        nut_.<span class="literal">min</span>(100);</span><br><span class="line">    &#125;</span><br><span class="line"></span><br><span class="line">    <span class="keyword">if</span> (debug)</span><br><span class="line">    &#123;</span><br><span class="line">        Info&lt;&lt; typeName &lt;&lt; ':' &lt;&lt; <span class="keyword">nl</span></span><br><span class="line">            &lt;&lt; <span class="string">"    max(Theta) = "</span> &lt;&lt; <span class="literal">max</span>(Theta_).value() &lt;&lt; <span class="keyword">nl</span></span><br><span class="line">            &lt;&lt; <span class="string">"    max(nut) = "</span> &lt;&lt; <span class="literal">max</span>(nut_).value() &lt;&lt; endl;</span><br><span class="line">    &#125;</span><br><span class="line">&#125;</span><br></pre></td></tr></table></figure></p>
<p>这里重点关注一下 <code>ThetaEqn</code> 的写法。<br><figure class="highlight css"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br></pre></td><td class="code"><pre><span class="line"><span class="tag">fvScalarMatrix</span> <span class="tag">ThetaEqn</span></span><br><span class="line">  (</span><br><span class="line">      1<span class="class">.5</span>*</span><br><span class="line">      (</span><br><span class="line">          <span class="rule"><span class="attribute">fvm</span>:<span class="value">:<span class="function">ddt</span>(alpha, rho, Theta_)</span><br><span class="line">        + fvm::<span class="function">div</span>(alphaRhoPhi, Theta_)</span><br><span class="line">        - fvc::<span class="function">Sp</span>(fvc::<span class="function">ddt</span>(alpha, rho) + fvc::<span class="function">div</span>(alphaRhoPhi), Theta_)</span><br><span class="line">      )</span><br><span class="line">    - fvm::<span class="function">laplacian</span>(kappa_, Theta_, <span class="string">"laplacian(kappa,Theta)"</span>)</span><br><span class="line">   ==</span><br><span class="line">      fvm::<span class="function">SuSp</span>(<span class="function">-</span>((PsCoeff*I) &amp;&amp; gradU), Theta_)</span><br><span class="line">    + (tau &amp;&amp; gradU)</span><br><span class="line">    + fvm::<span class="function">Sp</span>(-gammaCoeff, Theta_)</span><br><span class="line">    + fvm::<span class="function">Sp</span>(-J1, Theta_)</span><br><span class="line">    + fvm::<span class="function">Sp</span>(J2/(Theta_ + ThetaSmall), Theta_)</span><br><span class="line"> )</span></span>;</span><br></pre></td></tr></table></figure></p>
<p>对应的偏微分方程是，<br>$$<br>\frac{3}{2}\left [ \frac{\partial }{\partial t}(\varepsilon_s \rho_s \Theta)+ \nabla \cdot (\varepsilon_s \rho_s \Theta U_s) \right ] = (-P_s \mathrm{I} + \tau_s):\nabla U_s + \nabla \cdot (\kappa_s \nabla \Theta) - \gamma_s - J_s<br>$$<br>下面将公式与代码一一对应。代码里跟公式相比，多了一项<br>$$<br>-\Theta\frac{\partial }{\partial t}(\varepsilon_s \rho_s) - \Theta \nabla \cdot (\varepsilon_s \rho_s U_s)<br>$$<br>这样，相当于代码对应的前两项是<br>$$<br>\varepsilon_s \rho_s \frac{\partial \Theta}{\partial t} + \varepsilon_s \rho_s U_s \cdot \nabla \Theta<br>$$<br>这样做的目的仍不是很明确。</p>
<p>Laplacian 项不需多言，剩下的几项，全都当作了源项来处理。<br>$(-P_s \mathrm{I} + \tau_s)\:\nabla U_s$ 拆开成了两项，分别对应 <code>fvm::SuSp(-((PsCoeff*I) &amp;&amp; gradU), Theta_)</code> 和 <code>(tau &amp;&amp; gradU)</code> 。第一项，由于$P_s$是$\Theta$ 的函数，所以，在固相压力类中，返回的值是固相压力系数（固相压力除以颗粒温度），在这里将$\Theta$进行了隐式处理。而由于$\tau$与$\nabla U_s$与$\Theta$无关，所以就只当成一般的源项了。</p>
<p>$\gamma$项的处理与颗粒压力类似，也是只定义 <code>gammaCoeff</code> ，然后将$\Theta$作隐式处理。由于实际上$\gamma$的表达式<br>$$<br>\gamma = \frac{12(1-e^2)\varepsilon_s^2\rho g_0}{d_p \sqrt{\pi}} \Theta^{3/2}<br>$$<br>这里有关于$\Theta$的非线性项，程序里实际上是对$\Theta$作了部分隐式处理，从 <code>gammaCoeff</code>的定义可以看出来：<br>$$<br>\gamma_{Coeff} = \frac{12(1-e^2)\varepsilon_s^2\rho g_0}{d_p } \sqrt{\frac{\Theta}{\pi}}<br>$$</p>
<p>颗粒温度方程中的$J_s$ 也分成了两项来处理，$J_s = J_1 - J_2$<br>$J_1$的处理很简单，在公式里，$J_1 = 3\beta \Theta$，而在代码当中，变量 <code>J1</code> 定义为 <code>3*beta</code>，而 <code>Theta</code> 则作隐式处理。   </p>
<p>$J_2$也人为作了隐式处理，文献中$J_2$的表达式为<br>$$<br>J_2 = \frac{\beta^2 d_p |U_g-U_s|^2}{4\varepsilon_s\rho\sqrt{\pi \Theta}}<br>$$</p>
<p>程序里则将$J_2$除以$\Theta$以后作为系数，以实现对$\Theta$进行隐式处理。<br>注意这里要说明一下”Sp”和”SuSp”的区别，这个在”Programmer’s Guide”里有说明，此外这个<a href="https://openfoamwiki.net/index.php/HowTo_Adding_a_new_transport_equation" target="_blank" rel="external">网页</a>里也有说明。可是，为什么固相压力项要用”SuSp”，这个也暂时不明白。</p>
<p>计算颗粒温度的代数方程的具体公式这里不写了，可以参考文献。</p>
<p>最后，上面提到的由于颗粒摩擦作用产生的颗粒压力和颗粒粘度，只在颗粒体积分率很大的区域才需要启用，算例中需要设定一个值，对于颗粒相体积分率大于一个设定值时， <code>alphaMinFriction</code>，只有颗粒相体积分率超过这个值时，才会启用由于颗粒相摩擦而产生的压力和粘度。</p>

      
    </div>
    <footer class="article-footer">
      <a data-url="http://xiaopingqiu.github.io/2015/09/19/kineticTheoryModel/" data-id="cj2lay331002rkkmbrpyftd3a" class="article-share-link">Share</a>
      
      
  <ul class="article-tag-list"><li class="article-tag-list-item"><a class="article-tag-list-link" href="/tags/Code-Explained/">Code Explained</a></li><li class="article-tag-list-item"><a class="article-tag-list-link" href="/tags/OpenFOAM/">OpenFOAM</a></li></ul>

    </footer>
  </div>
  
    
<nav id="article-nav">
  
    <a href="/2015/09/19/kineticTheoryModelSubModels/" id="article-nav-newer" class="article-nav-link-wrap">
      <strong class="article-nav-caption">Newer</strong>
      <div class="article-nav-title">
        
          OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之 kineticTheoryModel 的子模型
        
      </div>
    </a>
  
  
    <a href="/2015/09/13/CPP-conversion-Of-derived-class-reference-to-base-class-type/" id="article-nav-older" class="article-nav-link-wrap">
      <strong class="article-nav-caption">Older</strong>
      <div class="article-nav-title">C++ 中派生类引用与基类引用的隐式转换</div>
    </a>
  
</nav>

  
</article>


<section id="comments">
  <!-- 多说评论框 start -->
  <div class="ds-thread" data-thread-key="post-kineticTheoryModel" data-title="OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之 kineticTheoryModel" data-url="http://xiaopingqiu.github.io/2015/09/19/kineticTheoryModel/"></div>
  <!-- 多说评论框 end -->
  <!-- 多说公共JS代码 start (一个网页只需插入一次) -->
  <script type="text/javascript">
  var duoshuoQuery = {short_name:"xiaopingqiu"};
	(function() {
		var ds = document.createElement('script');
		ds.type = 'text/javascript';ds.async = true;
		ds.src = (document.location.protocol == 'https:' ? 'https:' : 'http:') + '//static.duoshuo.com/embed.js';
		ds.charset = 'UTF-8';
		(document.getElementsByTagName('head')[0] 
		 || document.getElementsByTagName('body')[0]).appendChild(ds);
	})();
	</script>
  <!-- 多说公共JS代码 end -->
</section>

</section>
        
          <aside id="sidebar">
  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Categories</h3>
    <div class="widget">
      <ul class="category-list"><li class="category-list-item"><a class="category-list-link" href="/categories/C/">C++</a><span class="category-list-count">2</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/DEM/">DEM</a><span class="category-list-count">1</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/OpenFOAM/">OpenFOAM</a><span class="category-list-count">44</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/Paraview/">Paraview</a><span class="category-list-count">5</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/swak4Foam/">swak4Foam</a><span class="category-list-count">1</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/test/">test</a><span class="category-list-count">2</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/vim/">vim</a><span class="category-list-count">1</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Tags</h3>
    <div class="widget">
      <ul class="tag-list"><li class="tag-list-item"><a class="tag-list-link" href="/tags/Boundary-conditions/">Boundary conditions</a><span class="tag-list-count">6</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/C/">C++</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/CentOS/">CentOS</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Code-Explained/">Code Explained</a><span class="tag-list-count">29</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/LES/">LES</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/LIGGGHTS/">LIGGGHTS</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/ODE/">ODE</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/OpenFOAM/">OpenFOAM</a><span class="tag-list-count">20</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Postprocessing/">Postprocessing</a><span class="tag-list-count">9</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Preprocessing/">Preprocessing</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/RTS/">RTS</a><span class="tag-list-count">3</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/TIL/">TIL</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Windows/">Windows</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/fvOptions/">fvOptions</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/groovyBC/">groovyBC</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/paraview/">paraview</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/test/">test</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/thermophysicalModels/">thermophysicalModels</a><span class="tag-list-count">5</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/turbulence-model/">turbulence model</a><span class="tag-list-count">7</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/vim/">vim</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/wall-functions/">wall functions</a><span class="tag-list-count">4</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Tag Cloud</h3>
    <div class="widget tagcloud">
      <a href="/tags/Boundary-conditions/" style="font-size: 15.56px;">Boundary conditions</a><a href="/tags/C/" style="font-size: 11.11px;">C++</a><a href="/tags/CentOS/" style="font-size: 10px;">CentOS</a><a href="/tags/Code-Explained/" style="font-size: 20px;">Code Explained</a><a href="/tags/LES/" style="font-size: 10px;">LES</a><a href="/tags/LIGGGHTS/" style="font-size: 10px;">LIGGGHTS</a><a href="/tags/ODE/" style="font-size: 10px;">ODE</a><a href="/tags/OpenFOAM/" style="font-size: 18.89px;">OpenFOAM</a><a href="/tags/Postprocessing/" style="font-size: 17.78px;">Postprocessing</a><a href="/tags/Preprocessing/" style="font-size: 11.11px;">Preprocessing</a><a href="/tags/RTS/" style="font-size: 12.22px;">RTS</a><a href="/tags/TIL/" style="font-size: 10px;">TIL</a><a href="/tags/Windows/" style="font-size: 10px;">Windows</a><a href="/tags/fvOptions/" style="font-size: 11.11px;">fvOptions</a><a href="/tags/groovyBC/" style="font-size: 10px;">groovyBC</a><a href="/tags/paraview/" style="font-size: 10px;">paraview</a><a href="/tags/test/" style="font-size: 11.11px;">test</a><a href="/tags/thermophysicalModels/" style="font-size: 14.44px;">thermophysicalModels</a><a href="/tags/turbulence-model/" style="font-size: 16.67px;">turbulence model</a><a href="/tags/vim/" style="font-size: 10px;">vim</a><a href="/tags/wall-functions/" style="font-size: 13.33px;">wall functions</a>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Archives</h3>
    <div class="widget">
      <ul class="archive-list"><li class="archive-list-item"><a class="archive-list-link" href="/archives/2017/05/">五月 2017</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/08/">八月 2016</a><span class="archive-list-count">8</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/06/">六月 2016</a><span class="archive-list-count">5</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/05/">五月 2016</a><span class="archive-list-count">3</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/04/">四月 2016</a><span class="archive-list-count">12</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/03/">三月 2016</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/12/">十二月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/11/">十一月 2015</a><span class="archive-list-count">2</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/10/">十月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/09/">九月 2015</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/08/">八月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/06/">六月 2015</a><span class="archive-list-count">2</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/05/">五月 2015</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/04/">四月 2015</a><span class="archive-list-count">2</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Recents</h3>
    <div class="widget">
      <ul>
        
          <li>
            <a href="/2017/05/12/duoshuoAnnouncement/">多说评论系统将停止提供服务</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewScritps/">Paraview 脚本一例</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewCamera/">Paraview 中有关 Camera 的操作两例</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewCustomFilter/">Paraview 中创建 Custom Filter</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewStreamLineOnSlice/">在 Paraview 中画截面上的流线</a>
          </li>
        
      </ul>
    </div>
  </div>

  
</aside>
        
      </div>
      <footer id="footer">
  
  <div class="outer">
    <div id="footer-info" class="inner">
      &copy; 2017 Giskard Q.<br>
      Powered by <a href="http://hexo.io/" target="_blank">Hexo</a>
    </div>
  </div>
</footer>

<script src="//dn-lbstatics.qbox.me/lbservice/busuanzi/2.0/busuanzi.mini.js"/></script>

    </div>
    <nav id="mobile-nav">
  
    <a href="/" class="mobile-nav-link">Home</a>
  
    <a href="/archives" class="mobile-nav-link">Archives</a>
  
    <a href="/atom.xml" class="mobile-nav-link">Rss</a>
  
    <a href="/about" class="mobile-nav-link">About</a>
  
</nav>
    

<script src="//ajax.googleapis.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>


  <link rel="stylesheet" href="/fancybox/jquery.fancybox.css" type="text/css">
  <script src="/fancybox/jquery.fancybox.pack.js" type="text/javascript"></script>


<script src="/js/script.js" type="text/javascript"></script>

<script>
  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
  })(window,document,'script','//www.google-analytics.com/analytics.js','ga');

  ga('create', 'UA-62501686-1', 'auto');
  ga('send', 'pageview');

</script>

  </div>
<!-- mathjax config similar to math.stackexchange -->

<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    tex2jax: {
      inlineMath: [ ['$','$'], ["\\(","\\)"] ],
      processEscapes: true
    }
  });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {
        skipTags: ['script', 'noscript', 'style', 'textarea', 'pre', 'code']
      }
    });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Queue(function() {
        var all = MathJax.Hub.getAllJax(), i;
        for(i=0; i < all.length; i += 1) {
            all[i].SourceElement().parentNode.className += ' has-jax';
        }
    });
</script>

<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
</body>
</html>